Warnings to the user regarding the execution of the cells will appear in the red boxes:
How?
This notebooks uses Script of Scripts (SoS) to exchange data between Octave and Python subkernels. This is achieved by the %get magic command placed at the very beginning of the Python3 cells. For example:
%get var_in_octave --from Octave
Following this %get magic command, var_in_octave variable will be also available to Python3 cells.
.c implementations for calculating density compensation to use them in Octave%use Octave
disp('Executed automatically');
disp('Cell to load octave packages. To show, select the cell and click arrow up icon in the toolbar.');
pkg load image
pkg load optim
# Load python packages
print('Executed automatically');
print('Cell to import python modules. To show, select the cell and click arrow up icon in the toolbar.');
import plotly.plotly as py
import plotly.graph_objs as go
import plotly_express as px
from plotly import tools
import numpy as np
import ipywidgets as widgets
import math
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=True)
# This is a Python3 cell to create an interactive figure.
# Here we use %get magic function to migrate variables from the Octave workspace into Python
print('Executed automatically');
print('Cell to define some Plotly functions. To show, select the cell and click arrow up icon in the toolbar.');
def heatmap_trace(inp, name, xlen, ylen, clrmax, clrmin, clrscale):
trace = go.Heatmap( z =inp,
y = list(range(xlen-1)),
x = list(range(ylen-1)),
colorscale=clrscale,
showscale = True,
zmax=clrmax,
zmin=clrmin,
colorbar=dict(
tickfont=dict(
size=14,
color='white'
)),
name = name);
return trace
def heatmap_trace2(inp, name, xlen, ylen, clrscale):
trace = go.Heatmap( z =inp,
y = list(range(xlen-1)),
x = list(range(ylen-1)),
colorscale=clrscale,
showscale = False,
name = name);
return trace
% Mex c files for gridding by Brian Hargreaves and Philip Beatty
% http://mrsrl.stanford.edu/~brian/gridding/
disp('Executed automatically');
disp('Cell to MEX c functions for dcf. To show, select the cell and click arrow up icon in the toolbar.');
mex gridlut_mex.c
mex calcdcflut_mex.c
load('/tmp/rrsg_challenge/brain_radial_96proj_12ch.mat');
whos % Show variables in the current scope
rawdata = permute(rawdata,[4,3,2,1]);
trajectory = permute(trajectory,[3,2,1]);
[~,nFE,nSpokes,nCh] = size(rawdata);
whos
%get rawdata --from Octave
%get nFE --from Octave
%get nSpokes --from Octave
%get nCh --from Octave
print('Click on this cell and execute to pass variables from Octave to Python')
# This is a data exchange cell Octave --> Python3
# This cell plots rawdata from each channel
# The execution of this cell takes a bit longer than conventional 2D plots (e.g. matplotlib) as each datapoint is represented on a heatmap for interactivity using Plotly.
rawdata = np.squeeze(rawdata)
fig = tools.make_subplots(rows=2, cols=6, print_grid=False, horizontal_spacing = 0.02, vertical_spacing = 0.02)
traces = []
iter = 0
for ii in range(2):
for zz in range(6):
cur_trace = heatmap_trace(np.log(1+abs(rawdata[:,:,iter])), 'Channel: '+ str(iter+1), nSpokes, nFE, 0.0001, -0, 'Viridis')
fig.append_trace(cur_trace, ii+1, zz+1)
iter += 1
fig['layout'].update(height=600, width=800, title=dict(text='<b>Raw data from each channel</b> <br> <i>Hover to see channel number and data points.</i>',font=dict(color="#ffc107")),paper_bgcolor='#000000')
for ii in range(12):
exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')
exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')
# plot(fig,filename='rawdata.html')
print('Execute this cell to regenerate the following figure.')
iplot(fig)
from IPython.display import IFrame
IFrame(src='./rawdata.html', width=820, height=600,scrolling='no')
# First we will parse data on Octave
# Use bart to obtain root sum of squares of the rawdata over channels
rd = real(rawdata) + 1i*imag(rawdata);
rd = bart('rss 8', rd);
clr = squeeze(log(rd)); clear rd;
trajx = squeeze(trajectory(1,:,:));
trajy = squeeze(trajectory(2,:,:));
disp('Execute this cell before the follwing (if needed)');
%get trajx --from Octave
%get trajy --from Octave
%get clr --from Octave
# The reason trajx, trajy and trajz are not vectorized in the Octave cell above is that,
# currently SoS data exchange gets stuck while transferring large vectors.
trajx = np.reshape(np.array(trajx),(nFE*nSpokes))
trajy = np.reshape(np.array(trajy),(nFE*nSpokes))
clr = np.reshape(np.array(clr),(nFE*nSpokes))
trac = go.Scatter(
x = trajx,
y = trajy,
mode='markers',
marker=dict(
color = clr, #set color equal to a variable
colorscale= 'Viridis',
showscale=True,
colorbar=dict(
tickfont=dict(
size=14,
color='white'
))
),
)
layout = go.Layout(
autosize=False,
width=800,
height=800,
paper_bgcolor='#000000',
plot_bgcolor='#000000',
xaxis = go.layout.XAxis(
gridcolor= '#283442',
tickfont=dict(
size=14,
color='white'
)
),
yaxis = go.layout.YAxis(
gridcolor= '#283442',
tickfont=dict(
size=14,
color='white'
)
),
hovermode = 'closest',
title=dict(text='<b>Raw data projected on k-space trajectory</b> <br> <i>Hover to see the locations of each sample</i>',font=dict(color="#ffc107"))
)
fig = go.Figure(data=[trac], layout=layout)
#plot(fig, filename='kspace.html')
print('Execute this cell (after previous one) to reproduce the following figure.')
iplot(fig)
from IPython.display import IFrame
IFrame(src='./kspace.html', width=820, height=820,scrolling='no')
% Perform NUFFT to interpolate data onto cartesian grid.
% -d denotes dimension (x:y:z, which is 300X300X1)
% -i sets the transform type to inverse
% -l enables L2 regularization
% -t enables Toeplitz embedding
% trajectory is the k-space locations of the acquired samples
% rawdata is the sampled raw data
im = bart('nufft -d300:300:1 -i -l -t',trajectory,rawdata);
% Transform sensitivity maps to cartesian k-space using FFT
% -u denotes unitary transform
% 7 sets the bitmask level
im_ks = bart('fft -u 7', im);
% For details regarding BART's ECALIB, please see the implementation notes section.
calib = bart('ecalib -m1 -I',im_ks);
%get calib --from Octave
print('Transfer coil sensitivities to Python.')
# This is a data exchange cell. Octave --> Python3
calib = np.squeeze(calib)
fig = tools.make_subplots(rows=3, cols=4, print_grid=False, horizontal_spacing = 0.02, vertical_spacing = 0.02)
traces = []
iter = 0
for ii in range(3):
for zz in range(4):
mx = np.max(np.log(1+abs(calib[:,:,iter])))
mn = np.min(np.log(1+abs(calib[:,:,iter])))
cur_trace = heatmap_trace(np.log(1+abs(calib[:,:,iter])), 'Channel: '+ str(iter+1), 300, 300, mn,mx , 'Viridis')
fig.append_trace(cur_trace, ii+1, zz+1)
iter += 1
fig['layout'].update(height=750, width=800, title=dict(text='<b>Sensitivity profile of each channel</b> <br> <i>Hover to see channel number and data points.</i>',font=dict(color="#ffc107")),paper_bgcolor='#000000')
for ii in range(12):
exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels=False,showline = False, ticks = \'\')')
exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showline=False, showticklabels = False, ticks = \'\')')
# plot(fig,filename='calibrations.html')
print('Execute to renegerate coil sensitivity figure.');
iplot(fig)
from IPython.display import IFrame
IFrame(src='./calibrations.html', width=820, height=820,scrolling='no')
% Use gridded data for SENSE recon --> BART PICS
bart_SENSE = bart('pics -l2', im_ks, calib);
% Use non-cartesian SENSE recon --> BART PICS
bart_SENSE2 = bart('pics -t',trajectory, rawdata, calib);
%get bart_SENSE --from Octave
%get bart_SENSE2 --from Octave
fig = tools.make_subplots(rows=1, cols=2, print_grid=False,vertical_spacing = 0.02)
fig.append_trace(heatmap_trace2(abs(bart_SENSE), 'BART pics -l2', 300, 300,'Greys'),1,1)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2), 'BART pics -t', 300, 300,'Greys'),1,2)
fig['layout'].update(height=400, width=600, title=dict(text='<b>Cartesian vs non-cartesian SENSE in BART</b>',font=dict(color="#ffc107")),paper_bgcolor='#000000')
for ii in range(2):
exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')
exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showticklabels = False, ticks = \'\')')
#plot(fig,filename='bart_cart_non.html')
print('Execute this cell to reproduce following figure.')
iplot(fig)
from IPython.display import IFrame
IFrame(src='./bart_cart_non.html', width=800, height=420,scrolling='no')
% Please see HelperFunctions folder for subSample.m
[outRD_x2, outTR_x2] = subSample(rawdata,trajectory,2,nSpokes);
[outRD_x3, outTR_x3] = subSample(rawdata,trajectory,3,nSpokes);
[outRD_x4, outTR_x4] = subSample(rawdata,trajectory,4,nSpokes);
% BART non-cartesian sense outputs
bart_SENSE2_x2 = bart('pics -t',outTR_x2, outRD_x2, calib);
bart_SENSE2_x3 = bart('pics -t',outTR_x3, outRD_x3, calib);
bart_SENSE2_x4 = bart('pics -t',outTR_x4, outRD_x4, calib);
%get bart_SENSE --from Octave
%get bart_SENSE2_x2 --from Octave
%get bart_SENSE2_x3 --from Octave
%get bart_SENSE2_x4 --from Octave
print('Execut to pass variables from Octave to Pyhton.')
# Data exchange cell. Octave --> Python3
fig = tools.make_subplots(rows=2, cols=2, print_grid=False,vertical_spacing = 0.02, horizontal_spacing=0.04)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2), 'Original', 300, 300,'Greys'),1,1)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2_x2), '2X', 300, 300,'Greys'),1,2)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2_x3), '3X', 300, 300,'Greys'),2,1)
fig.append_trace(heatmap_trace2(abs(bart_SENSE2_x4), '4X', 300, 300,'Greys'),2,2)
fig['layout'].update(height=600, width=600, title=dict(text='<b>Comparison of BART non-cartesian SENSE outputs <br> at different undersampling rates</b>',font=dict(color="#ffc107")),paper_bgcolor='#000000')
for ii in range(4):
exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, showline=False, zeroline = False, showticklabels = False, ticks = \'\')')
exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showline=False, showticklabels = False, ticks = \'\')')
#plot(fig,filename='sub_bart_compare.html')
print('Execute to reproduce the following figure');
iplot(fig)
from IPython.display import IFrame
IFrame(src='./sub_bart_compare.html', width=620, height=620,scrolling='no')
Reminder: Figure-1 from the paper
Implementation of iterative image reconstruction.
Conjugate gradient (CG) iteration is controlled by the central CG process. It is initialized by MR data originating from N receiver channels (1,2,3,…N), acquired with an arbitrary k‐space trajectory.
Separately for each channel, these data undergo processing similar to conventional gridding reconstruction, i.e., sampling density correction (D), and resampling along a Cartesian grid, followed by FFT (FT1). The resulting images are individually multiplied by complex conjugate coil sensitivity and summed. After subsequent intensity correction (I), the sum image represents the vector a as defined by Eq. [25]. After initialization with a, the CG process iteratively calculates a progression of images, which converges towards exact reconstruction. For each iteration step, a current residuum image vector needs to be multiplied by the matrix IEH DEI. This is performed by the loop which starts from the CG box. After initial intensity correction (I), the processing is continued separately for each receiver coil. First, the intensity corrected residuum image is multiplied by individual coil sensitivity (Si). The results are transformed into k‐space by FFT and resampled along the experimental k‐space trajectory (FT2), resulting in a set of multiple‐coil k‐space data similar to that obtained experimentally. The following steps are equivalent to those carried out with the original data, yielding an intermediate image, which is fed back into the CG process. Here a refined approximate solution is calculated, which serves for further refinement by continued iteration.
As soon as the current approximation is sufficiently accurate, it is output and undergoes final intensity correction and k‐space filtering.
Implementation notes
This section relates the present implementation with the Figure-1 using the same notations.
FT1
operation is implemented using BART's NUFFT:This function transforms non-cartesian k-space data (rawdata and trajectory) into image domain (im):
im = bart('nufft -d x:y:z -i -l -t',trajectory,rawdata);
% -d denotes dimensions in x:y:z
% -i sets the transform type to inverse
% -l enables L2 regularization
% -t enables Toeplitz embedding
% trajectory is the k-space locations of the acquired samples
% rawdata is the sampled raw data
S$\gamma$
coil sensitivities are estimated using BART's ECALIB. See section 2.1This function creates sensitivity maps (stored in the calib variable) from the raw k-space data on cartesian grid (sens_maps_ks):
calib = bart('ecalib -m1 -I',sens_maps_ks);
% -I enables intensity correction
% -m1 sets number of maps to one
FT2
is implemented using BART's NUFFT:This function transforms image im back to a non-uniform k-space of (radial, which is defined by trajectory variable in this case) nu_ks:
nu_ks = bart('nufft',traj,im);
D
is the operation of scaling rawdata with density correction matrix.To obtain density correction factor dcf pertaining to the rawdata, submodules (written in C) provided by Brian Hargreaves's gridding functions are used:
dcf = calcdcflut(trajectory,300);
% trajectory is the k-space locations of the non-cartesian rawdata
% 300 is the recon matrix size (300x300)
Sum
is implemented using BART's RRS:This function simply performs sum of squares combination:
sum = bart('rss bitmask',multi_channel_input);
F
is the k-space filtering, defined in 3.6D
NOTE: This operation may take a few minutes.
A simpler method could have been applied here. Instead, I preferred using the submodules (written in C) provided by Brian Hargreaves's gridding functions to provide a good example of interoperability.
dcf = calcdcflut(trajectory,300);
dcf = reshape(dcf,[3 nFE nSpokes]);
I
# Root sum of square of sensitivities (estimated in section 2.1) from 12 channels
I = abs(bart('rss 8',calib.*conj(calib)));
%get I --from Octave
trace = heatmap_trace2(abs(I), 'I', 300, 300,'Viridis')
axis_template = dict(showgrid = False, zeroline = False,
linecolor = 'black', showticklabels = False,
ticks = '' )
layout = dict(height=500, width=500, title=dict(text='<b>Intensity correction matrix</b>',font=dict(color="skyblue")),paper_bgcolor='#000000',xaxis=axis_template,yaxis=axis_template)
figure = dict(data=[trace],layout=layout)
#plot(figure,filename='incor.html')
print('Execute to generate the following figure');
from IPython.display import IFrame
IFrame(src='./incor.html', width=500, height=500,scrolling='no')
E
function E = opE(inp,S,traj,I)
Input arguments
Output
The following cell executes automatically to ensure that function definition exists.
function E = opE(inp,S,traj,I)
inp = inp.*I; % Scale with the intensity correction matrix, as operation I precedes operation E (see Fig. 1)
tmp = S.*inp; % Multiply the intensity corrected image with coil sensitivities. This will produce one image per channel stored in tmp variable
E = bart('nufft',traj,tmp); % Transform back to the non-uniform k-space (see implementation notes for details)
end
EH
function EH = opEH(dcf,inp,S,traj,I)
Input arguments
Output
The following cell executes automatically to ensure that function definition exists.
function EH = opEH(dcf,inp,S,traj,I)
tmp = bart('nufft -d300:300:1 -i -l -t',traj,inp.*sqrt(dcf(1,:,:))); % NUFFT to image domain (see implementation notes)
Sstar = conj(S(:,:,1,:)); % Complex conjugate of sensitivity profiles (see Fig. 1)
tmp2 = bart('rss 8 ',tmp.*Sstar); % Images scaled by complex conjugate of sensitivity profiles and SOS combined
EH = tmp2.*I; % Scale output image by intensity correction
end
F
The following cell executes automatically to ensure that function definition exists.
function out = opF(im,beta,r,N,I)
im = im./I;
%imagesc(im);
im(isnan(im)==1) = 0;
f_k = zeros(N,N);
f_k(N/2+1,N/2+1) = 1;
f_k = bwdist(f_k);
f_k = 0.5 + 1/pi.*atan(beta.*((r-abs(f_k))/r));
f_k = fftshift(f_k);
im_k = fft2(im);
out = abs(ifft2(ifftshift(im_k.*f_k)));
end
The following cell executes automatically to ensure that function definition exists.
function [b,deltas] = cg_solve(a,I,S,dcf,maxstep,trajectory)
p = a;
r = a;
b = zeros(300,300);
deltas = zeros(maxstep,1);
for ii=1:maxstep
disp(['Iteration -->' num2str(ii)]);
delta = r(:)'*r(:)/(a(:)'*a(:));
disp(delta);
deltas(ii) = delta;
E = opE(p,S,trajectory);
q = opEH(dcf,E,S,trajectory,I);
% dot(r,r) is equivalent to r(:)'*r(:). Used dot for easy reading.
term = dot(r,r)/dot(p,q);
b = b + term*p;
rprev = r;
r = r - term*q;
term2 = dot(r,r)/dot(rprev,rprev);
p = r + term2*p;
end
end
The following cell executes automatically to ensure that function definition exists.
function [im,deltas] = main_sense(rawdata, calib, trajectory,I,dcf,maxiter)
a = opEH(dcf,rawdata,calib,trajectory,I);
[im,deltas] = cg_solve(a,I,calib,dcf,maxiter,trajectory);
im = opF(im,100,40,300,I);
end
[im,deltas] = main_sense(rawdata, calib, trajectory,I,dcf,10);
[im_x2,deltas_x2] = main_sense(outRD_x2, calib, outTR_x2,I,dcf(:,:,1:2:end),10);
[im_x3,deltas_x3] = main_sense(outRD_x3, calib, outTR_x3,I,dcf(:,:,1:3:end),10);
[im_x4,deltas_x4] = main_sense(outRD_x4, calib, outTR_x4,I,dcf(:,:,1:4:end),10);
%get deltas --from Octave
%get deltas_x2 --from Octave
%get deltas_x3 --from Octave
%get deltas_x4 --from Octave
iterations = [1,2,3,4,5,6,7,8,9,10];
sub2_1 = go.Scatter(
x = iterations,
y = np.log10(deltas),
name = 'R=1',
text = 'R=1',
hoverinfo = 'x+y+text',
line=dict(color='rgb(0,114,189)'),
)
sub2_2 = go.Scatter(
x = iterations,
y = np.log10(deltas_x2),
name = 'R=2',
text = 'R=2',
hoverinfo = 'x+y+text',
line=dict(color='rgb(217,83,25)'),
)
sub2_3 = go.Scatter(
x = iterations,
y = np.log10(deltas_x3),
name = 'R=3',
text = 'R=3',
hoverinfo = 'x+y+text',
line=dict(color='rgb(237,177,32)'),
)
sub2_4 = go.Scatter(
x = iterations,
y = np.log10(deltas_x4),
name = 'R=4',
text = 'R=4',
hoverinfo = 'x+y+text',
line=dict(color='rgb(126,47,142)'),
)
fig = tools.make_subplots(rows=1, cols=2, print_grid=False)
fig.append_trace(sub2_1, 1, 1)
fig.append_trace(sub2_2, 1, 1)
fig.append_trace(sub2_3, 1, 1)
fig.append_trace(sub2_4, 1, 1)
iplot(fig)
%get im --from Octave
%get im_x2 --from Octave
%get im_x3 --from Octave
%get im_x4 --from Octave
fig = tools.make_subplots(rows=2, cols=2, print_grid=False,vertical_spacing = 0.02, horizontal_spacing=0.04)
fig.append_trace(heatmap_trace2(abs(im), 'Original', 300, 300,'Greys'),1,1)
fig.append_trace(heatmap_trace2(abs(im_x2), '2X', 300, 300,'Greys'),1,2)
fig.append_trace(heatmap_trace2(abs(im_x3), '3X', 300, 300,'Greys'),2,1)
fig.append_trace(heatmap_trace2(abs(im_x4), '4X', 300, 300,'Greys'),2,2)
fig['layout'].update(height=600, width=600, title=dict(text='<b>Comparison of iterative algorithm <br> at different undersampling rates</b>',font=dict(color="#ffc107")),paper_bgcolor='#000000')
for ii in range(4):
exec('fig[\'layout\'][\'xaxis' + str(ii+1) + '\'].update(showgrid = False, showline=False, zeroline = False, showticklabels = False, ticks = \'\')')
exec('fig[\'layout\'][\'yaxis' + str(ii+1) + '\'].update(showgrid = False, zeroline = False, showline=False, showticklabels = False, ticks = \'\')')
#plot(fig,filename='sub_bart_compare.html')
iplot(fig)